home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 December / macformat-019.iso / Shareware City / Developers / Gnuplot 3.5 for Macintosh / demo / airfoil.dem < prev    next >
Encoding:
Text File  |  1993-11-03  |  4.7 KB  |  172 lines  |  [TEXT/MSWD]
  1. #
  2. # $Id: airfoil.dem%v 3.38.2.45 1993/01/11 04:09:41 woo Exp woo $
  3. #
  4. # This demo shows how to use bezier splines to define NACA four
  5. # series airfoils and complex variables to define Joukowski
  6. # Airfoils.  It will be expanded after overplotting in implemented
  7. # to plot Coefficient of Pressure as well.
  8. #        Alex Woo, Dec. 1992
  9. #
  10. # The definitions below follows: "Bezier presentation of airfoils",
  11. # by Wolfgang Boehm, Computer Aided Geometric Design 4 (1987) pp 17-22.
  12. #
  13. #                Gershon Elber, Nov. 1992
  14. #
  15. # m = percent camber
  16. # p = percent chord with maximum camber
  17. pause 0  "NACA four series airfoils by bezier splines"
  18. pause 0  "Will add pressure distribution later with Overplotting"
  19. mm = 0.6
  20. # NACA6xxx
  21. thick = 0.09  
  22. # nine percent  NACAxx09
  23. pp = 0.4
  24. # NACAx4xx
  25. # Combined this implies NACA6409 airfoil
  26. #
  27. # Airfoil thickness function.
  28. #
  29. set xlabel "NACA6409 -- 9% thick, 40% max camber, 6% camber"
  30. x0 = 0.0
  31. y0 = 0.0
  32. x1 = 0.0
  33. y1 = 0.18556
  34. x2 = 0.03571
  35. y2 = 0.34863
  36. x3 = 0.10714
  37. y3 = 0.48919
  38. x4 = 0.21429 
  39. y4 = 0.58214
  40. x5 = 0.35714
  41. y5 = 0.55724
  42. x6 = 0.53571
  43. y6 = 0.44992
  44. x7 = 0.75000
  45. y7 = 0.30281
  46. x8 = 1.00000
  47. y8 = 0.01050
  48. #
  49. # Directly defining the order 8 Bezier basis function for a faster evaluation.
  50. #
  51. bez_d4_i0(x) =     (1 - x)**4
  52. bez_d4_i1(x) = 4 * (1 - x)**3 * x
  53. bez_d4_i2(x) = 6 * (1 - x)**2 * x**2
  54. bez_d4_i3(x) = 4 * (1 - x)**1 * x**3
  55. bez_d4_i4(x) =                  x**4
  56.  
  57. bez_d8_i0(x) =      (1 - x)**8
  58. bez_d8_i1(x) =  8 * (1 - x)**7 * x
  59. bez_d8_i2(x) = 28 * (1 - x)**6 * x**2
  60. bez_d8_i3(x) = 56 * (1 - x)**5 * x**3
  61. bez_d8_i4(x) = 70 * (1 - x)**4 * x**4
  62. bez_d8_i5(x) = 56 * (1 - x)**3 * x**5
  63. bez_d8_i6(x) = 28 * (1 - x)**2 * x**6
  64. bez_d8_i7(x) =  8 * (1 - x)    * x**7
  65. bez_d8_i8(x) =                   x**8
  66.  
  67.  
  68. m0 = 0.0
  69. m1 = 0.1
  70. m2 = 0.1
  71. m3 = 0.1
  72. m4 = 0.0
  73. mean_y(t) = m0 * mm * bez_d4_i0(t) + \
  74.         m1 * mm * bez_d4_i1(t) + \
  75.         m2 * mm * bez_d4_i2(t) + \
  76.         m3 * mm * bez_d4_i3(t) + \
  77.         m4 * mm * bez_d4_i4(t)
  78.  
  79. p0 = 0.0
  80. p1 = pp / 2
  81. p2 = pp
  82. p3 = (pp + 1) / 2
  83. p4 = 1.0
  84. mean_x(t) = p0 * bez_d4_i0(t) + \
  85.         p1 * bez_d4_i1(t) + \
  86.         p2 * bez_d4_i2(t) + \
  87.         p3 * bez_d4_i3(t) + \
  88.         p4 * bez_d4_i4(t)
  89.  
  90. z_x(x) = x0 * bez_d8_i0(x) + x1 * bez_d8_i1(x) + x2 * bez_d8_i2(x) + \
  91.      x3 * bez_d8_i3(x) + x4 * bez_d8_i4(x) + x5 * bez_d8_i5(x) + \
  92.      x6 * bez_d8_i6(x) + x7 * bez_d8_i7(x) + x8 * bez_d8_i8(x)
  93.  
  94. z_y(x, tk) = \
  95.    y0 * tk * bez_d8_i0(x) + y1 * tk * bez_d8_i1(x) + y2 * tk * bez_d8_i2(x) + \
  96.    y3 * tk * bez_d8_i3(x) + y4 * tk * bez_d8_i4(x) + y5 * tk * bez_d8_i5(x) + \
  97.    y6 * tk * bez_d8_i6(x) + y7 * tk * bez_d8_i7(x) + y8 * tk * bez_d8_i8(x)
  98.  
  99. #
  100. # Given t value between zero and one, the airfoild curve is defined as
  101. # 
  102. #            c(t) = mean(t1(t)) +/- z(t2(t)) n(t1(t)),
  103. #
  104. # where n is the unit normal to the mean line. See the above paper for more.
  105. #
  106. # Unfortunately, the parametrization of c(t) is not the same for mean(t1)
  107. # and z(t2). The mean line (and its normal) can assume linear function t1 = t,
  108. #                                                     -1
  109. # but the thickness z_y is, in fact, a function of z_x  (t). Since it is
  110. # not obvious how to compute this inverse function analytically, we instead
  111. # replace t in c(t) equation above by z_x(t) to get:
  112. # 
  113. #            c(z_x(t)) = mean(z_x(t)) +/- z(t) n(z_x(t)),
  114. #
  115. # and compute and display this instead. Note we also ignore n(t) and assumes
  116. # n(t) is constant in the y direction,
  117. #
  118.  
  119. airfoil_y1(t, thick) = mean_y(z_x(t)) + z_y(t, thick)
  120. airfoil_y2(t, thick) = mean_y(z_x(t)) - z_y(t, thick)
  121. airfoil_y(t) = mean_y(z_x(t))
  122. airfoil_x(t) = mean_x(z_x(t))
  123. set nogrid
  124. set nozero
  125. set parametric
  126. set xrange [-0.1:1.1]
  127. set yrange [-0.1:.7]
  128. set trange [ 0.0:1.0]
  129. set title "NACA6409 Airfoil"
  130. plot airfoil_x(t), airfoil_y(t) title "mean line" w l 2, \
  131.      airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l 1, \
  132.      airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l 1
  133. pause -1 "Press Return"
  134. mm = 0.0
  135. pp = .5
  136. thick = .12
  137. set title "NACA0012 Airfoil"
  138. set xlabel "12% thick, no camber -- classical test case"
  139. plot airfoil_x(t), airfoil_y(t) title "mean line" w l 2, \
  140.      airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l 1, \
  141.      airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l 1
  142. pause -1 "Press Return"
  143. set title ""
  144. set xlab ""
  145. set key
  146. set parametric
  147. set samples 100
  148. set isosamples 10
  149. set data style lines
  150. set function style lines
  151. pause 0  "Joukowski Airfoil using Complex Variables" 
  152. set title "Joukowski Airfoil using Complex Variables" 0,0
  153. set time
  154. set yrange [-.2 : 1.8]
  155. set trange [0: 2*pi]
  156. set xrange [-.6:.6]
  157. zeta(t) = -eps + (a+eps)*exp(t*{0,1})
  158. eta(t) = zeta(t) + a*a/zeta(t)
  159. eps = 0.06
  160. a =.250
  161. set xlabel "eps = 0.06 real"
  162. plot real(eta(t)),imag(eta(t))
  163. pause -1 "Press Return"
  164. eps = 0.06*{1,-1}
  165. set xlabel "eps = 0.06 + i0.06"
  166. plot real(eta(t)),imag(eta(t))
  167. pause -1 "Press Return"
  168. set title ""
  169. set xlabel ""
  170.  
  171.  
  172.